##################
###Reduced Form###
##################

####Baseline
#Subset na omit data
ttt <- na.omit(dis.dat[,c("outbreak_event","chicken_meat_tonnes1000000","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000")])
#Estimate baseline model
fs1.rf = lm(outbreak_event ~ comb.np.use_tonnes1000000+t+factor(ccode), data = dis.dat)
fs1.rf.sum <- coeftest(fs1.rf, vcov = vcovCL(fs1.rf, cluster = ~ccode))

####Full
#Subset na omit data
ttt2 <- na.omit(dis.dat[,c("outbreak_event","outbreak_event","chicken_meat_tonnes1000000","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])

##Run full model
fs2.rf = lm(outbreak_event ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = dis.dat)
fs2.rf.sum <- coeftest(fs2.rf, vcov = vcovCL(fs2.rf, cluster = ~ccode))

#####Zero only
##Create a combined meat indicator to assess zero-to-irelevant food production
dis.dat$all_meat_tonnes1000000 <- dis.dat$cattle_meat_tonnes1000000+dis.dat$chicken_meat_tonnes1000000+dis.dat$pigs_meat_tonnes1000000
summary(dis.dat$all_meat_tonnes1000000)
##Subset little-to-irrelevant food production countries
zero.only <- subset(dis.dat, all_meat_tonnes1000000<0.050)
#Remove NAs
zero.only$logGDPpc[is.na(zero.only$logGDPpc)] <- 0
zero.only$comb.np.use_tonnes1000000[is.na(zero.only$comb.np.use_tonnes1000000)] <- 0
#Keep only relevant varialbles 
subset.dat2 <- na.omit(zero.only[,c("outbreak_event","flus_all","all_meat_tonnes1000000","cattle_meat_tonnes1000000","chicken_meat_tonnes1000000","pigs_meat_tonnes1000000","comb.np.use_tonnes1000000","t","lagoutbreak_event", "lagflus_all","logpopulation","logGDPpc","ccode", "year")]) 
#Estimate model
test.iv.rf = lm(outbreak_event ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+
                  logpopulation+logGDPpc+factor(ccode), data = subset.dat2)
summary(test.iv.rf)
test.iv.rf.sum <- coeftest(test.iv.rf, vcov = vcovCL(test.iv.rf, cluster = ~ccode))
test.iv.rf.sum

#####Rest of the sample
##Subset little-to-irrelevant food production countries
non.zeros <- subset(dis.dat, all_meat_tonnes1000000>0.050)
#Keep only relevant varialbles 
subset.dat3 <- na.omit(non.zeros[,c("outbreak_event","flus_all","all_meat_tonnes1000000","cattle_meat_tonnes1000000","chicken_meat_tonnes1000000","pigs_meat_tonnes1000000","comb.np.use_tonnes1000000","t","lagoutbreak_event", "lagflus_all","logpopulation","logGDPpc","ccode", "year")]) 
#Estimate model
test.iv.rf.2 = lm(outbreak_event ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+
                  logpopulation+logGDPpc+factor(ccode), data = subset.dat3)
summary(test.iv.rf.2)
test.iv.rf.2.sum <- coeftest(test.iv.rf.2, vcov = vcovCL(test.iv.rf.2, cluster = ~ccode))
test.iv.rf.2.sum


#####################################################
#####################################################
###Minimal and Regular/high Production First Stage###
#####################################################
#####################################################

###None production sample
test.iv.di1 = lm(all_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+
                   logpopulation+logGDPpc+factor(ccode), data = subset.dat2 )
summary(test.iv.di1)
test.iv.di.sum1 <- coeftest(test.iv.di1, vcov = vcovCL(test.iv.di1, cluster = ~ccode))
test.iv.di.sum1

##Rest of sample
test.iv.di2 = lm(all_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+
                   logpopulation+logGDPpc+factor(ccode), data = subset.dat3)
summary(test.iv.di2)
test.iv.di.sum2 <- coeftest(test.iv.di2, vcov = vcovCL(test.iv.di2, cluster = ~ccode))
test.iv.di.sum2


##Export OLS table to Latex
stargazer(test.iv.di.sum1,test.iv.di.sum2)
#Obs and diagnostics
stargazer(fs1.rf,fs2.rf, test.iv.rf, test.iv.rf.2,test.iv.di1,test.iv.di2)



#################
#################
###First Stage###
#################
#################

##############
###Baseline###
##############
#Subset NA data
ttt <- na.omit(dis.dat[,c("outbreak_event","cattle_meat_tonnes1000000","chicken_meat_tonnes1000000","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000")])
#Beef
fs1.cattle = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+factor(ccode), data = ttt)
#Chicken meat
fs1.chicken = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+factor(ccode), data = ttt)
#Pork
fs1.pigs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+factor(ccode), data = ttt)

##########
###Full###
##########
##Subset NAs
ttt2 <- na.omit(dis.dat[,c("outbreak_event","cattle_meat_tonnes1000000","chicken_meat_tonnes1000000","pigs_meat_tonnes1000000","t","ccode","comb.np.use_tonnes1000000","lagoutbreak_event","logpopulation","logGDPpc")])

##Beef
fs2.cattle = lm(cattle_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt2)
##Chicken meat
fs2.chicken = lm(chicken_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt2)
##Pork
fs2.pigs = lm(pigs_meat_tonnes1000000 ~ comb.np.use_tonnes1000000+t+lagoutbreak_event+logpopulation+logGDPpc+factor(ccode), data = ttt2)

##Export OLS table to Latex
stargazer(fs1.cattle, fs1.chicken, fs1.pigs, fs2.cattle, fs2.chicken, fs2.pigs)




